실습: 경사 충격파 관계식#
강좌: 수치해석 프로젝트
이론#
수직충격파#
초음속 유동이 물체이 부딛히거나 유동의 방향이 꺽이는 경우 충격파 또는 팽창파가 발생한다. 수직 방향으로 충격파가 발생하는 경우 충격파 전/후 마하수는 다음과 같은 Prandtl 관계를 만족한다.
여기서 \(M^*\) 는 특성 마하수로 Sonic point일 때 음속으로 마하수 \(M\) 과는 아래 관계를 만족한다.
여기서 \(\gamma\) 는 비열비로 일반 공기는 1.4이다.
이 식을 이용하면 Prandtl 관계식은 다음과 같이 정리할 수 있다.
충격파 전/후 물성치 변화는 다음과 같다.
충격파 전/후에는 전압력 손실이 발생한다. 충격파를 제외한 영역은 등엔트로피 유동으로 가정할 수 있다. 즉 마하수가 \(M\) 인 유동의 정압 (Static Pressure, \(p\)) 대비 전압력 (Total Pressure, \(p_0\)) 는 다음 관계를 갖는다.
전압력 손실은 충격파 전/후의 전압력을 각각 구한 후 그 감소 비율이다.
경사충격파#
쐐기 각도 \(\theta\) 인 경우 충격파의 파각은 \(\beta\)로 발생한다. 이 때 충격파에 수직인 성분은 수직충격파 관계식을 따르고, 수평 방향 속도는 같다. 즉
Fig. 3 Shock Wave#
\(\theta-\beta-M\) 에 대해서는 아래 관계식을 만족한다.
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
plt.style.use('ggplot')
plt.rcParams['figure.dpi'] = 150
def tbm(beta, M, gamma=1.4):
"""
Theta-Beta-M relation
Parameters
----------
beta : float
Angle of wave (rad)
M : float
Mach number
gamma : float
specific heat ratio
Returns
-------
theta : float
Angle of wedge (rad)
"""
# tangent(theta)
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
# tan^{-1}(theta)
return np.arctan(tangent)
M = 1.1
beta = np.deg2rad(10)
theta = np.rad2deg(tbm(beta, M))
print('Theta={} @ M={}, beta={}'.format(theta, M, beta))
beta = np.deg2rad(70)
theta = np.rad2deg(tbm(beta, M))
print('Theta={} @ M={}, beta={}'.format(theta, M, beta))
Theta=-66.15222865219691 @ M=1.1, beta=0.17453292519943295
Theta=1.0317298907919357 @ M=1.1, beta=1.2217304763960306
from scipy import optimize
def tbm_curve(M, gamma=1.4, n=31):
# Construct equation for beta
f = lambda x : tbm(x, M, gamma)
# Compute root of tbm w.r.t beta
sol = optimize.root_scalar(f, bracket=[0, np.pi/2])
beta_min = sol.root
# Divide betas on (beta_min, 90)
beta = np.linspace(beta_min, np.pi/2, n)
# Compute thetas
theta = tbm(beta, M)
return theta, beta
# Plot for Mach numbers
for M in [1.1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.5, 3, 4, 6, 15]:
# Get theta-beta curve (in degree)
theta, beta = np.rad2deg(tbm_curve(M))
# Plot each curve
plt.plot(theta, beta, color='gray')
# Get maximum theta and beta
idx = theta.argmax()
theta_max = theta[idx]
beta_loc = beta[idx]
# Add text
plt.text(theta_max, beta_loc, "M={}".format(M), fontsize='x-small')
# Add labels using Latex symbols
plt.xlabel(r"$\theta$")
plt.ylabel(r"$\beta$")
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
/tmp/ipykernel_71941/1007713943.py:20: RuntimeWarning: divide by zero encountered in double_scalars
tangent = 2/np.tan(beta) * (M**2*np.sin(beta)**2 - 1) / (M**2*(gamma+np.cos(2*beta))+2)
Text(0, 0.5, '$\\beta$')

실습#
수직충격파#
마하수에 따른 전압력을 계산하는 코드를 작성하시오.
충격파 전 마하수 M1 일 때 충격파 후 마하수 및 밀도 및 압력비를 구하는 코드를 작성하시오
def p0_p(M, gamma=1.4):
#...
return ...
def normal_shock(M1, gamma=1.4):
# ...
return M2, rho2, p2, p0_loss
경사충격파#
충격파 전 마하수 M1, 쐐기 각도가 theta 일 때 \(\theta-\beta-M\) 비선형 방정식을 해석하시오.
\(f(\beta, M) \rightarrow \tan(\theta\)) 를 구하는 코드를 만드시오.
\(\beta\)가 달라질 때 \(-\tan(\theta)\) 를 최소화하는 조건과 그 때 \(\theta\)를 구하시오.
\(\beta\) 를 미지수 x 라 하고, 주어진 \(M, \theta\) 에 대해서 \(f(x, M) - \tan(\theta)=0\) 인 방정식의 해를 수치적으로 구하시오.
충격파 후 마하수, 밀도, 압력비 및 파각을 구하는 코드를 작성하시오.
Hint) 파각에 수직방향 마하수 \(M_{n1}, M_{n2}\)
\[ M_{n1} = \frac{M_1}{\sin(\beta)}, M_{n2} = \frac{M_2}{\sin(\beta-\theta)} \]
def tangent_theta(beta, M, gamma=1.4):
# Calculuate theta-beta-M relation
return t_theta
def theta_max(M, gamma=1.4):
# function to compute -tangent_theta(x, M, gamma=1.4) for given M
# minimize the function with very small x0
return theta
def beta_week(M, theta, gamma=1.4):
# make equation f(x,M) - tan(theta) =0
# root finding
return beta
def oblique_shock(M1, theta, gamma=1.4):
# Solve theta-beta-M
# Compute Mn1
# Solve normal shock for Mn1
# Calculate M2
return M2, rho2, p2, p0_loss beta